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ABSTRACT 

This  study  presents  a  methodology  for  computing 
stochastic  sensitivities  with  respect  to  the  design  variables, 
which  are  the  mean  values  of  the  input  correlated  random 
variables.  Assuming  that  an  accurate  surrogate  model  is 
available,  the  proposed  method  calculates  the  component 
reliability,  system  reliability,  or  statistical  moments  and  their 
sensitivities  by  applying  Monte  Carlo  simulation  (MCS)  to  the 
accurate  surrogate  model.  Since  the  surrogate  model  is  used, 
the  computational  cost  for  the  stochastic  sensitivity  analysis  is 
negligible.  The  copula  is  used  to  model  the  joint  distribution 
of  the  correlated  input  random  variables,  and  the  score 
function  is  used  to  derive  the  stochastic  sensitivities  of 
reliability  or  statistical  moments  for  the  correlated  random 
variables.  An  important  merit  of  the  proposed  method  is  that  it 
does  not  require  the  gradients  of  performance  functions,  which 
are  known  to  be  erroneous  when  obtained  from  the  surrogate 
model,  or  the  transformation  from  X-space  to  U-space  for 
reliability  analysis.  Since  no  transformation  is  required  and 
the  reliability  or  statistical  moment  is  calculated  in  X-space, 
there  is  no  approximation  or  restriction  in  calculating  the 
sensitivities  of  the  reliability  or  statistical  moment.  Numerical 
results  indicate  that  the  proposed  method  can  estimate  the 
sensitivities  of  the  reliability  or  statistical  moments  very 
accurately,  even  when  the  input  random  variables  are 
correlated. 

KEYWORDS 


Stochastic  Sensitivity  Analyses,  Score  Functions,  Monte  Carlo 
Simulation,  Copula,  Surrogate  Model,  RBDO 

1.  INTRODUCTION 

Reliability-based  design  optimization  (RBDO)  and 
reliability-based  robust  design  optimization  (RBRDO)  have 
been  widely  applied  to  various  engineering  applications  such 
as  stamping  [1,2],  vehicle  design  with  durability  [3,4],  and 
noise,  vibration,  and  harshness  (NVF1)  analysis  [5,6],  where 
accurate  sensitivities  of  performance  functions  are  available.  If 
accurate  sensitivities  are  available  in  a  complex  physical 
system,  then  the  most  probable  point  (MPP)-based  reliability 
analysis,  which  includes  the  First  Order  Reliability  Method 
(FORM)  [7-10],  the  Second  Order  Reliability  Method 
(SORM)  [11,12],  and  the  MPP -based  Dimension  Reduction 
Method  (DRM)  [13-15],  can  be  used  for  approximately 
assessing  the  reliability  of  the  system,  which  is  used  as  a 
probabilistic  constraint  in  both  RBDO  and  RBRDO. 
Furthermore,  the  first  and  second  statistical  moments  of  the 
performance  function,  which  are  used  in  the  objective  function 
of  RBRDO,  are  approximated  using  the  first-order  Taylor 
series  expansion  [16-17]  or  the  mean-based  DRM  [18-20]. 

However,  for  engineering  applications  where  accurate 
sensitivities  of  performance  functions  are  not  available,  such 
as  drivetrain,  crashworthiness,  micro  or  nano  mechanics,  or 
fluid-structure  interaction,  the  MPP -based  reliability  analysis, 
which  uses  the  sensitivities  of  performance  functions  to  find 
the  MPP,  cannot  be  directly  used.  Instead,  surrogate  models 


have  been  widely  used  to  carry  out  design  optimization  for  the 
engineering  applications  where  sensitivities  are  unavailable 
[21-23].  Once  an  accurate  surrogate  model  is  available  for  the 
design  optimization,  the  direct  Monte  Carlo  simulation  (MCS) 
[24]  can  be  applied  to  the  accurate  surrogate  model  to  estimate 
the  reliability  or  statistical  moments  of  the  system  with 
negligible  computational  burden.  The  reliability  or  statistical 
moment  will  be  called  the  probabilistic  response  in  this  paper. 
To  use  the  probabilistic  response  obtained  using  the  MCS  for 
the  design  optimization,  sensitivities  of  the  probabilistic 
response  are  still  required,  which  can  be  obtained  using  the 
finite  difference  method  (FDM)  or  using  the  sensitivities  of 
the  surrogate  model.  However,  since  the  probabilistic  response 
is  obtained  from  the  MCS,  the  FDM  may  require  a  significant 
number  of  samples  to  obtain  accurate  sensitivities  due  to  the 
statistical  noise  of  the  MCS.  On  the  other  hand,  even  if  the 
surrogate  model  is  very  accurate  for  the  response  value,  the 
sensitivities  obtained  from  the  surrogate  model  are  known  to 
be  inaccurate,  and  accordingly  it  is  not  a  good  idea  to  use  them 
for  the  design  optimization. 

The  main  objective  of  this  paper  is  to  propose  a  new 
sampling-based  RBDO  for  correlated  random  inputs  [25-27], 
which  does  not  require  obtaining  the  sensitivities  of  the 
performance  functions  and  their  sensitivities  from  the 
surrogate  model.  Instead,  stochastic  sensitivity  analysis  using 
the  score  function,  which  was  proposed  for  the  independent 
random  variables  [28]  or  correlated  Gaussian  random 
variables  [29],  is  used  to  derive  the  sensitivities  of  the 
probabilistic  response  with  correlated  random  variables.  The 
proposed  sensitivity  analysis  does  not  require  the 
transformation  from  the  original  design  space  to  the  standard 
normal  space,  which  means  that  there  is  no  approximation  or 
restriction  in  calculating  the  sensitivities  of  the  probabilistic 
response.  Since  the  generation  of  the  accurate  surrogate  model 
is  beyond  the  scope  of  this  paper,  it  will  not  be  explained  here. 
For  this,  the  dynamic  Kriging  (D-Kriging)  method  developed 
in  Ref.  31  is  used  in  this  paper.  Numerical  examples 
demonstrate  that  the  sensitivities  developed  in  the  paper  agree 
very  well  with  the  sensitivities  obtained  from  the  FDM  with 
250  million  samples  for  the  MCS. 

The  paper  is  organized  into  four  main  parts.  The  first  part. 
Section  2,  explains  how  to  mathematically  express  the 
component,  the  system  probability  of  failure,  and  the 
statistical  moments  and  their  sensitivities.  The  second  part. 
Section  3,  shows  how  to  derive  the  score  function  for 
statistically  independent  or  correlated  random  variables.  The 
third  part,  Section  4,  shows  the  formulation  of  the  sampling- 
based  RBDO  and  RBRDO.  The  last  part.  Section  5, 
demonstrates  with  numerical  examples  the  accuracy  of  the 
proposed  sensitivities  compared  with  the  FDM  results.  Finally, 
the  proposed  method  is  integrated  with  the  D-Kriging  method 
to  carry  out  RBDO  and  compared  with  four  other  RBDO 
methods  to  validate  and  point  out  the  merit  of  the  proposed 
method. 

2.  RELIABILITY,  STATISTICAL  MOMENTS  AND 
SENSITIVITY 

2.1  Reliability  and  Statistical  Moments 

A  reliability  analysis,  for  both  the  component  and  the 
system  level,  involves  calculation  of  the  probability  of  failure. 


denoted  by  PF,  which  is  defined  using  a  multi-dimensional 
integral 

PAv)  -  FIX  e  nF]  =  jRV 4f  (x)fx(x;y)dx  =  2?[/n,  (X)]  (1) 

where  v|/  is  a  vector  of  distribution  parameters,  which  usually 
includes  the  mean  (p)  and/or  standard  deviation  (a)  of  the 
random  input  X  =  ,  /*[•]  represents  a 

probability  measure,  Q,  is  the  failure  set,  /X(x;i|/)  is  a  joint 
probability  density  function  (PDF)  of  X,  and  is[«]  represents 
the  expectation  operator.  The  failure  set  is  defined  as 
Q,  =  |  x  :  G  ( x )  >  0}  for  component  reliability  analysis  of  the 

ith  constraint  function  G;(x),  and  ClF  =  {x:U;>,«>  oj  and 

={x:n,=iG-(x)>o}  for  series  system  and  parallel 

system  reliability  analysis  of  NC  performance  functions, 
respectively  [29,30].  /Qf  ( x)  in  Eq.  (1)  is  called  an  indicator 

function  and  defined  as 

fl,  xeQr 

/n  (x)=  F  (2) 

F  [0,  otherwise 

In  this  paper,  since  the  mean  of  X,  p  =  { }T ,  is 

used  as  a  design  vector,  the  vector  of  distribution  parameters 
v|/  is  simply  replaced  with  p  for  the  derivation  of  the 
sensitivity. 

In  a  fashion  similar  to  Eq.  (1),  the  qth  statistical  moment 
of  a  performance  function  H(x)  is  defined  as 

m,(p)  =  E[ff<l(X)\  =  j\  H«(x)fx(x;n)dx,  (3) 

and  thus,  Eqs.  (1)  and  (3)  can  be  written  in  a  generalized  form 
as 

/7(n)  ee  £[g(X)]  =  j\,  g(x)/x(x;n)7/x  (4) 

which  is  called  a  probabilistic  response  [29].  In  Eq.  (4),  /?(p) 
and  g(x)  represent  PF{yi)  and  /,,  (x)  ,  respectively,  for 
reliability  analysis,  and  h(\i)  and  g(x)  represent  and 

Hr,(x) ,  respectively,  for  statistical  moment  analysis. 

2.2  Stochastic  Sensitivity  Analysis 

The  sensitivity  of  the  probabilistic  response  hip,)  with 
respect  to  //,  is  considered  in  this  section.  For  the  derivation 
of  the  sensitivity,  the  following  four  assumptions,  which  are 
known  as  the  regularity  conditions,  are  required  [28,29]. 

1.  The  joint  PDF  /X(x;p)  is  continuous. 

2.  The  mean  jut  e  M;  c  R,  i  =  where  M;  is  an  open 

interval  on  R . 


3. 


4. 


The  partial  derivative 


g/xC^lO 

8 Hi 


exists  and  is  finite  for  all 


x  and  //, .  In  addition,  /?(p)  is  a  differentiable  function  of 


p. 

There  exists  a  Lebesgue  integrable  dominating  function 
r(x)  such  that 


g(x) 


g/x(x;p) 

dju, 


<r(x) 


(5) 


for  all  p. 

With  the  four  assumptions  satisfied,  taking  the  partial 
derivative  of  Eq.  (4)  with  respect  to  fi,  yields 

g/i(n)  _  _g_  r  g(x)  f  (\-  p)t/x  (6) 

8/u,  dju,  JrK 


and  the  differential  and  integral  operators  can  be  interchanged 
due  to  Assumption  4  and  the  Lebesgue  dominated 
convergence  theorem  [29,32],  giving 


dh{]i) 

dju, 


-\,sWdM ail* 

gln/x(x;p) 
dju, 

,  gin  /X(x;n) 


=  Lv  £(x)  /x (x;  p)Jx 

g(x)' 


=  £ 


g//, 


(7) 


since  g(x)  is  not  a  function  of  ju,  .  The  partial  derivative  of  the 
log  function  of  the  joint  PDF  in  Eq.  (7)  with  respect  to  ji,  is 
known  as  the  first-order  score  function  for  u,  and  is  denoted 
as 


^(x;n); 


gln/x(x;p) 

dju, 


(8) 


Therefore,  it  is  required  to  know  the  first-order  score 
function  to  derive  the  sensitivity  of  the  probabilistic  response, 
which  is  either  the  reliability  or  the  statistical  moments.  The 
derivation  of  the  first-order  score  function  for  independent  and 
correlated  random  variables  will  be  shown  in  the  subsequent 
section. 


3.  SENSITIVITY  ANALYSIS  BY  SCORE  FUNCTION 
3.1  For  Independent  Random  Variables 

Consider  a  random  input  X=  whose 

components  are  statistically  independent  random  variables. 
Then,  the  joint  PDF  of  X  is  expressed  as  multiplication  of  its 
marginal  PDFs  as 


fx  (x;  f)  =  n,=i  fx,  ;a) 


(9) 


where  fx  ( x, ;  ju, )  is  the  marginal  PDF  corresponding  to  the  i,h 

random  variable  X;.  Therefore,  for  statistically  independent 
random  variables,  the  first-order  score  function  for  u,  is 
expressed  as 


(10) 

'  o//,  5//,. 

Since  the  marginal  PDF  and  cumulative  distribution 
function  (CDF)  are  available  analytically  as  listed  in  Table  1, 
where  <t>(.)  and  </>(•)  are  the  standard  normal  CDF  and  PDF, 
respectively,  given  by 

®(W)  =  £_  M&dZ  =  -2=  £  exp  y%,  (11) 

the  derivation  of  Eq.  (10)  is  straightforward  for  normal, 
lognormal,  and  Gumbel  distributions. 


Table  1.  Marginal  PDF,  CDF,  and  Its  Parameters 


PDF,  fx(x) 

CDF,  Fx(x) 

Parameters 

Normal 

1  -0.5[^f 

c  a 

jU,  (7 

V2/r<7 

{  <r  J 

Log- 

nonnal 

1  0.5| 

c  a 

hvj 

a2  =ln[l  +  (— )2] , 

\[27TX& 

F 

ji  =  ln(//)  -0.5a2 

Gumbel 

-a(x-v) 

ae 

0.577  n 

^  =  V+  ’C7=  rz 

Weibull 

V  V 

,  -o* 

l-e  v 

//  =  vr(i+i), 

k 

v2[r(i +y)~  r2(i  +  i)] 
k  k 

Uniform 

1 

- ,  a  <  x  <b 

b  —  a 

x- a 

b  —  a 

a  +  b  b -  a 

jU  = - ,  (7  =  . — 

2  x/l2 

However,  for  the  case  of  Weibull  distribution,  the 
derivation  is  not  straightforward  since  two  distribution 
parameters,  k  and  v,  are  coupled  as  shown  in  Table  1.  The 
score  function  for  Weibull  distribution  is  written  as 


■s®(x;iO  =  — 

dii, 


In  —  +  (k,  —1)  In  —  —  (—)*' 
v,  v v, 


1  dk,  1  dv,  dk,  ,  x, 

= - - '-  +  — Lln  — 

k,  djuj  v,  dju,  dju,  v , 

(k,  - 1)  8  v,  ( x ,  f  8k,  ^  k,  8  v, " 

v,  dJui  v,  ^  8/u,  v,  v,  dju, 


(12) 


which  requires  the  calculation  of  — —  and  — L  since  k,  and  v,- 

d/j,  dju, 

are  functions  of  Li .  As  shown  in  Table  1,  both  ^-L-  and  — L 

dju,  dju, 

require  the  evaluation  of  the  gamma  function  f  defined  as 


/•oo  1  / 

r(z)  =  J  t  e  dt  and  its  derivative.  By  taking  the  partial 

derivative  on  two  parameters  of  Weibull  distribution  and  using 

dk. 

the  derivative  of  the  gamma  function,  we  can  obtain  — and 

djUi 


dv 

djUf 


as 


where  X  and  |i  are  vectors  of  the  random  variables  and  its 
means,  respectively,  without  the  i,h  component.  Since  two 
integral  limits,  a  and  b,  are  functions  of  the  differential 
variable  /q,  the  differential  and  integral  operators  cannot  be 
interchanged  directly.  Instead,  the  Leibniz  integral  rule  [33] 
gives  a  formula  for  differentiation  of  a  definite  integral  whose 
limits  are  functions  of  the  differential  variable,  which  is  given 
by 


dVi 

o Mi 

dkt 

d/ui 

'.  +  i' 

V  Kj 


]/ .  /•  co  , 

(- [  tk‘e~'\ntdt 

k2  Jo 


2vX 


1  +  - 


2Vj 

~kT 


2 

,  k,  —t 


t  ‘  e  1  In  tdt 


i 

2  ^ 


(13) 


and  inserting  Eq.  (13)  into  Eq.  (12)  yields  the  first-order  score 
function  for  Weibull  distribution.  Table  2  summarizes  the 

dju. 

first-order  score  functions  for  four  marginal  PDFs.  — L  and 

d/Ji 

-^L  for  lognormal  distribution  in  Table  2  can  be  easily 

dfUi 

derived  using  the  definition  of  Jlj  and  a,  shown  in  Table  1, 

and  and  ^-L  for  Weibull  distribution  in  Table  2  are 
dju.  d jui 

shown  in  Eq.  (13). 


Table  2.  First-Order  Score  Function  for  u  for  Independent 
Random  Variables 


Marginal 

Distribution 


First-Order  Score  Function,  s[|)(x;|i) 


Nonnal 


X,  - 


Log-normal 

Gumbel 

Weibull 


1  da  1  [  In  x,  -  u 

— ‘-  +  —^  - 

cr  8Jui  a~  ^  a: 


a.  -  a.e 


-  SMi  n  -  ^ 

ct(— L  +  Onx,-^)— l 

OjUl  o/ui 


1  dki  1  dvt  +  dki  {ki  - 1)  dvt 

kj  dUi  d/ni  djni  v. 


V- 


dk,  k,  dv, 

8jui  vt  vj  d/ui 


v,.  dfi, 


Finally,  the  uniform  distribution,  as  shown  in  Table  1,  is 
not  continuous  on  K ,  which  makes  the  joint  PDF  in  Eq.  (9) 
discontinuous.  This  discontinuity  of  the  uniform  distribution 
violates  the  first  assumption  explained  in  Section  2.2.  Thus, 
the  score  function  cannot  be  used  if  a  random  variable  follows 
uniform  distribution. 

Assuming  X;  follows  the  uniform  distribution  and  all  the 
components  of  X  are  statistically  independent  random 
variables,  then  the  sensitivity  of  the  probabilistic  response 
with  respect  to  //,  shown  in  Eq.  (6)  can  be  written  as 


dft(n) 

dju, 


—  f  f 

dju JrAM 


*(«) 

a(ft) 


g(x)fx  (*,• ;  Mi )./v  (x;  j i)<Mx  ( 14) 


d_ 

dz 


fl'>f(x,z)dx  =  jhi')^dx  +  f(b(z),z)^- 

Ja(z)  J  a(z)  rj7 


db 


°(z)  dz 


dz 


-f(a(z),z)^- 

dz 


(15) 


Applying  the  Leibniz  integral  rule  to  Eq.  (14)  yields 


dhj lQ 

6 Mi 


g(x,b)-g(x,a) 

b-a 


/x(x;ii  )dx 


g(X,fc)-g(X,q) 

b-a 


(16) 


which  is  the  sensitivity  of  the  probabilistic  response  with 
respect  to  /q  when  X;  follows  uniform  distribution. 

3.2  For  Correlated  Input  Random  Variables 

Consider  a  bivariate  correlated  random  input 

X  =  jx  , ALj1  .  Then,  the  joint  PDF  of  X  is  expressed  as 
[26,27] 


f  ,  ,3  C(u,v;  G) 

/x  (x;  p)  =  — — - fx,  (x, ;  m,  )fx,  (Xj  ,Mj) 

dudv  1  (17) 

=  Cm  (; u ,  v;  9)fx  (x. ;  //,.  )fXj  (x , ;  jUj ) 

where  C  is  a  copula  function,  u  =  Fx  (xi;Juj)  and 

v=Fx  (x j ;  //  )  are  CDFs  for  X;  and  Xj,  respectively,  and  9  is 

the  correlation  coefficient.  Table  3  lists  commonly  used 
copula  functions  [26,27]. 

_ Table  3.  Commonly  Used  Copula  Functions _ 


Copula 

Type 


Copula  Function,  C ( u , \\9) 


Clayton 
AMH 

Frank 
FGM 

Gaussian  £'WJ 
Independent 


(«-*+»-*• 


if 


uv 


1  — 0(l  — w)(l  — v) 

fin[l+(e-“-l)(e-*-l)/(e-“-l)] 
uv  +  Guv(\-u)(\-v ) 


r®~‘M 


2^VT? 


rexp 


f  9 

2 Osw  —  s2  -  w2 


2(l-6>2) 


dsdw 


uv 


The  partial  derivative  of  the  copula  function  with  respect 
to  the  marginal  CDFs  u  and  v  is  called  the  copula  density  Frank 

function  and  denoted  as 


0 


eW+»)  _gW+v)  +g»("+v) 


+  1 


du 

5/1, 


c(u,v  |  6)  = 


d2C(u,v;0) 

dudv 


=  C  (m,v;6>) 


(18) 


and  displayed  in  Table  4.  As  shown  in  Table  4,  the  joint  PDF 
for  independent  random  variables  explained  in  Eq.  (9)  is  a 
special  case  of  Eq.  (17)  where  the  independent  copula  function 
is  used. 


Table  4.  Copula  Density  Functions 
Copula  Type  Copula  Density  Functions,  c(u,  v;  0) 


Clayton 

AMH 

Frank 

FGM 

Gaussian 

Independent 


(1  +  0)(uv)H1+m  (-1  +  ue  +  v-y2+e) 
1  +  02  (1  -  «)(l  -  v)  -  0(2  —  u  -  v-  uv ) 
[l-0(l-«)(l-v)]3 

0e(,<-'m)[ee  -l) 


1  +  <9(1  -  2m)(1  -  2v) 
i  |  »-'(«)■ 


‘V  2 

'  2<9®“1(w)®“1(v)  -  -  ®~‘(v)2 


exp 


2(1  -  02 ) 


Accordingly,  using  Eq.  (17),  the  first-order  score 
functions  in  Eq.  (8)  for  a  correlated  bivariate  input  are 
expressed  as 


FGM 


Gaussian 

Independent 


26(2v-l) 

l  +  0(l-2«)(l-2v) 


du 

djuj 


®  \u)  6®  ‘(v)-®  \u)  du 

^-\u))  +  <K®-\uW-02)\^, 


0 


Table  6.  Partial  Derivatives  of  Marginal  Distribution  with 
_ respect  to  //,■ _ 


Maiginal  Partial  Derivatives  of  Marginal  Distribution,  — — 

Distribution  5/1, 


Normal 


J*2.  1  if 

L~f^exp  “T 

V2ir<7,  2  ^ 


1  U-M, 


2 1  <j. 


r 


Log-nonnal  p  e  5ct,.  f  ln£-/i,  5/1,  |  (ln£-/i,)2  da, 


Gumbel 

Weibull 

Unifonn 


cr,  du, 

-a,  (Jf,  -V, 


djUi 


,h_  )*,• 

(  \ 

k‘  ( 

V,-  1 

x. 

\ 

dk ,  ^  x,  kj  dvt 

dju ,  v,  v,  y 

1 


b  —  a 


In  a  similar  way  to  that  explained  in  Section  3.1,  consider 
a  bivariate  correlated  random  input  X  =  j  A, ,  A,  J  where  X;  is 

assumed  to  follow  the  uniform  distribution.  Then,  the 
sensitivity  of  the  probabilistic  response  with  respect  to  /i,  for 
the  bivariate  correlated  random  input  can  be  written  as 


a), . ^  _  51n/x(x;fi)  _  d\nc{u,v,6)  i  gln/^lx,;/!,) 


(x;p)  = 


5/i, 


9ju, 


.(19) 


The  derivation  of  the  first  term  of  the  right-hand  side  of 
Eq.  (19)  is  straightforward  using  Table  4  and  listed  in  Table  5, 
and  the  second  term  of  the  right-hand  side  of  Eq.  (19)  is 
identical  to  Eq.  (10),  so  it  can  be  obtained  from  Table  2.  One 
can  see  from  Table  5  that  Eq.  (19)  is  identical  to  Eq.  (10)  if  the 
independent  copula  is  used.  In  Table  5,  the  partial  derivative 

Du 

of  the  marginal  CDF  with  respect  to  /i„  -  is  easily  obtained 

0/1, 

from  Table  1  and  is  shown  in  Table  6. 


Table  5.  Log-derivative  of  Copula  Density  Function 


Copula 

Type 

51nc(«,v;#) 

0/1, 

Clayton 

f  1  +  0  (20  +  1  )tT(1+9)> 
1  u  ire  +  v~e  - 1  j 

|  du 

1 5/i, 

AMH 

Q1  (l  -  v)  +  6(y  + 1) 

1 

m 

du 

1  -  0~  (1  -  «)(1  -  v)  —  0(2  —  u  -  v-  uv)  1 

-0(l-u)(l-v) 

0A 

&(x)c(m,  v;  6)  fx  (x, ;  /i,  )fXj  (x . ;  /i ,  )dxidxj 


Shjfi)  _ 

djui  dnt  JrJc,(a> 


(20) 


Applying  the  Leibniz  integral  rule  to  Eq.  (20)  yields 

dh(\i)  _  rbW  g(x j,xt)  dc(u,v,0) 

5/1,  b-a  0/i,  (21) 

g(X  ,b)c(X  b)-g(X  a)c(X a) 

+ - - - - - - - - - 1 

b-a 

and  Eq.  (21)  is  identical  with  Eq.  (16)  if  the  independent 
copula  is  used  since  c(u,  v;  0)  =  1 . 

4.  For  Both  Independent  and  Correlated  Input 
Random  Variables 

Consider  a  random  input  X  =  {X1,---,XWJT  where  M 

pairs  of  bivariate  correlated  random  variables  exist.  Then,  the 
joint  PDF  of  X  is  expressed  as 

/x  (x;  n)=n  r=i  cj  (u’ v;  )n  Zi  f*. ;  a  ) 


(22) 


from  Eqs.  (9)  and  (17).  Taking  the  partial  derivative  on  both 
sides  of  Eq.  (22)  yields  the  first-order  score  functions  for  a 
general  random  input  X  as 


(26) 


du 

d\nc(u,v;d)  51n  /v  (x,.;/v,.)  ~  ! 

=  - - - —  + - ! - . 

djUi  djuj 

Thus,  if  X;  is  a  correlated  random  variable,  Eq.  (23)  is 
identical  to  Eq.  (19),  and  if  X;  is  statistically  independent,  Eq. 
(23)  is  identical  to  Eq.  (10).  Similarly,  if  X;  follows  the 
uniform  distribution,  Eq.  (21)  can  be  used  to  calculate  the 
sensitivity  of  the  probabilistic  response  with  a  general  random 
input  X. 

5.  FORMULATION  OF  SAMPLING-BASED  RBDO 

The  mathematical  formulation  of  a  general  RBDO 
problem  is  expressed  as 

minimize  Cost(d) 

subject  to  P[G,.(X)  >  0]  <  Pj™ ,  i  =  \,-,NC  (24) 

dL  <  d  <  du,  d  e  R"'A'  and  X  e  R"n’ 

where  d  =  {<7,.  }T  =  p(X)  is  the  design  vector,  which  is  the 
mean  value  of  the  /V-dimensional  random  vector 
\={Xt ,  XN}J  ;  PjM  is  the  target  probability  of  failure 

for  the  Ith  constraint;  and  NC,  ndv,  and  nrv  are  the  number  of 
probabilistic  constraints,  design  variables,  and  random 
variables,  respectively.  The  mathematical  formulation  of  a 
general  RBRDO  is  given  by 

minimize  /  (/jh  ,  ) 

subject  to  P[G, (X)  >  0]  <  P'* ,  i  =  l,  —  ,NC  (25) 
dL  <  d  <  du,  de  R"*’  and  X  e  R",v 

where  fiH  and  a2H  are  the  mean  value  and  variance  of  the 
performance  function  H(X),  respectively. 

To  carry  out  RBDO  and  RBRDO  using  Eqs.  (24)  and 
(25),  respectively,  it  is  required  to  know  the  function  value 
and  its  sensitivities  at  a  given  design.  However,  in  most 
engineering  applications,  it  is  very  difficult,  if  not  impossible, 
to  obtain  accurate  sensitivities.  For  engineering  applications 
where  accurate  sensitivities  are  not  available,  suiTogate 
models  have  been  widely  used  to  carry  out  design 
optimization.  Once  an  accurate  surrogate  model  is  available 
for  the  design  optimizations,  the  MCS  can  be  applied  to 
estimate  the  reliability  or  statistical  moments  of  the  system 
with  negligible  computational  burden. 

Denote  the  surrogate  models  for  constraint  function  G;(X) 

and  performance  function  H(X)  as  Gj(X)  and  // ( X )  , 
respectively.  Then,  by  applying  the  MCS  to  the  surrogate 
model  G,(X) ,  the  probabilistic  constraints  in  Eqs.  (24)  and 
(25)  can  be  approximated  as 


PF  ,P[G,.(X)>0]  =  il4  (x(t))  <PTf 

x  k=i  F 


where  K  is  the  MCS  sample  size,  \<kl  is  the  k,h  realization  of 
X,  and  the  failure  set  Q,  for  the  surrogate  model  is  defined  as 

Q ,,  =  jx  :  G,  (x)  >  o|  .  Sensitivity  of  the  probabilistic  constraint 


in  Eqs.  (24)  and  (25)  is  obtained  using  the  score  function 
explained  in  Section  3  as 


where  ^(x  ,n)  is  obtained  using  Eq.  (23). 

In  a  similar  manner,  the  statistical  moments  in  Eq.  (25) 
can  be  approximated  as 

and  (28) 

A  k=x  A  k=x 

respectively,  and  their  sensitivities  are 

^-^^(x^.VV)  (29) 

8  jut  f  h 


|f4£"V*XV*>;|. .)-2^,  (30) 


respectively. 

As  shown  in  Eqs.  (27),  (29),  and  (30),  the  sensitivity 
calculation  using  the  score  function  and  MCS  does  not  require 
the  sensitivity  of  the  surrogate  models,  which  is  known  to  be 
inaccurate  even  if  the  surrogate  model  accurately 
approximates  the  function  values.  Furthermore,  the 
computation  of  the  sensitivity  using  the  score  function  does 
not  include  any  approximation  except  the  statistical  noise  due 
to  the  MCS,  which  can  be  avoided  using  a  sufficiently  large 
MCS  sample  size.  In  addition,  this  sensitivity  analysis  does 
not  require  the  transformation  from  the  original  design  space 
to  the  standard  normal  space,  which  usually  makes  the 
performance  function  become  highly  nonlinear,  especially 
when  the  random  input  follows  non-Gaussian  marginal 
distribution  and  is  correlated.  Therefore,  the  sensitivity 
analysis  using  the  score  function  and  MCS  will  be  very 
accurate  and  computationally  efficient  for  engineering 
applications  with  correlated  random  input  once  accurate 
suiTogate  models  are  available;  this  will  be  verified  through 
numerical  examples  in  the  subsequent  section. 

6.  NUMERICAL  EXAMPLES 

Numerical  studies  are  carried  out  in  this  section  to  verify 
the  stochastic  sensitivities  derived  in  Section  3  using  the  score 
function.  For  the  benchmark  sensitivity  to  test  the  accuracy  of 
the  proposed  method,  the  FDM  using  the  MCS  with  250 


million  samples  is  used.  The  stochastic  sensitivities  of  the 
component  probability  of  failure  and  statistical  moments  are 
compared  with  the  FDM  results  in  Section  5.1  and  Section  5.2, 
respectively.  Section  5.3  illustrates  how  the  proposed 
sensitivity  combined  with  an  accurate  surrogate  model  can  be 
used  to  solve  an  RBDO  problem. 

6.1  Sensitivities  of  Component  Probability  of  Failure 

Consider  a  2-D  highly  nonlinear  polynomial  function, 
which  was  studied  in  Ref.  15, 

G,  (X)  =  0.7361 +  (F-6)2  +(7-6)3  -0.6x(7-6)4  +Z  (31) 


where 


0.9063  0.42261  (X, 
0.4226  -0.9063 J[X2 


.  As  shown  in  Table  7, 


Xi  and  X2  follow  Weibull  and  normal  distributions, 
respectively.  The  Weibull  distribution  with  the  scale 
parameter  of  4  and  shape  parameter  of  15  has  a  mean  of 
3.8627  and  a  standard  deviation  of  0.3159.  Two  random 
variables  Xi  and  X2  are  assumed  to  be  correlated  with  each 
copula  shown  in  Table  3.  Except  for  the  independent  copula 
case  where  the  correlation  coefficient  is  always  zero,  the 
correlation  of  coefficient  of  x=0.3  is  used.  For  the  sensitivity 
calculation  of  the  component  probability  of  failure  using  the 
score  function,  the  sample  size  for  the  MCS  is  1  million, 
whereas  the  sensitivity  calculation  using  the  FDM  employs 
250  million  samples  for  the  MCS.  This  is  because  the  FDM 
uses  the  difference  of  two  probabilities  of  failure  obtained 
from  the  MCS,  which  could  include  relatively  large  statistical 
noise  if  a  small  number  of  samples  are  used. 


Table  7.  Properties  of  Random  Variables 


Random 

Variables 

Distribution 

Distribution  Parameters 

Xi 

Weibull 

v=4 

k=  15 

X2 

Normal 

iu=2.5 

cr=0.3 

Clayton 

9.7944 

9.7446 

AMH 

9.7148 

9.7354 

Frank 

8.9506 

8.9208 

FGM 

9.7167 

9.7421 

Gaussian 

7.4054 

7.4158 

Independent 

9.4483 

9.4296 

Table  9.  Sensitivity  of  Probability  of  Failure  w.r.t  //, 


Sensitivity 

Copula 

Score 

Function 

FDM  with 

perturbation  size  0.1% 

Clayton 

0.6485 

0.6513 

AMH 

0.6701 

0.6651 

Frank 

0.8149 

0.8321 

FGM 

0.9318 

0.9394 

Gaussian 

0.8347 

0.8165 

Independent 

0.0897 

0.0884 

To  verify  whether  the  proposed  stochastic  sensitivity 
analysis  works  for  high-dimensional  problems,  consider  a  4-D 
polynomial  function 


G,  (X)  =  145  - 1 1A,  +  X\  +  2  A,2  -  20A2 
-10A3+a32-21A4+2A42 


(32) 


The  properties  of  the  random  variables  in  Eq.  (32)  are 
shown  in  Table  10.  For  this  problem,  two  random  variables  X] 
and  X2  are  assumed  to  be  correlated  with  each  copula  shown 
in  Table  3,  and  two  random  variables  X3  and  X4  are  assumed 
to  be  statistically  independent.  The  relatively  large  correlation 
coefficient  (t=0.7)  for  X ,  and  X2  is  used  except  for  the 
independent  copula.  Flowever,  since  the  AMF1  copula  cannot 
deal  with  a  correlation  coefficient  larger  than  1/3  [27],  the 
AMF1  copula  is  excluded  in  this  example. 


Tables  8  and  9  show  the  comparison  of  sensitivities  with 
respect  to  the  mean  of  the  non-normal  random  variable  Xi  and 
the  normal  random  variable  X2,  respectively.  Both  tables 
illustrate  that  the  sensitivities  of  the  probability  of  failure 
obtained  using  the  score  function  in  Eq.  (27)  and  MCS  with  1 
million  samples  agree  very  well  with  the  sensitivities  obtained 
using  the  FDM  and  MCS  with  250  million  samples.  As 
mentioned  previously,  this  good  agreement  is  very  obvious 
since  the  sensitivity  derivation  explained  in  Section  3  does  not 
include  any  approximation  except  the  statistical  noise  due  to 
the  MCS.  Flence,  regardless  of  the  marginal  distribution  and 
copula  types  of  random  variables,  the  stochastic  sensitivity 
using  the  score  function  is  very  accurate,  and  it  is  also 
computationally  efficient  since  the  MCS  uses  a  relatively 
small  sample  size  compared  with  the  sensitivity  analysis  using 
the  FDM,  and  no  finite  difference  perturbation  size  is 
involved. 

Table  8.  Sensitivity  of  Probability  of  Failure  w.r.t  //, 

Sensitivity 

Copula  Score  FDM  with 

Function  perturbation  size  0.1% 


Table  10.  Properties  of  Random  Variables 


Random 

Variables 

Distribution 

Mean 

Standard 

Deviation 

Xi 

Normal 

5.0 

0.3 

x2 

Normal 

6.0 

0.3 

X3 

Normal 

7.0 

0.3 

X4 

Normal 

4.0 

0.3 

Table  11  compares  two  sensitivities  of  the  probability  of 
failure  with  respect  to  //, .  As  in  the  previous  example,  two 
sensitivities  agree  very  well  with  each  other.  Accordingly,  it 
can  be  said  that  the  proposed  stochastic  sensitivity  analysis  for 
the  probability  of  failure  is  very  accurate,  regardless  of  the 
type  of  performance  functions  and  copulas. 


Table  11.  Sensitivity  of  Probability  of  Failure  w.r.t  //, 


Sensitivity 

Copula 

Score 

FDM  with 

Function 

perturbation  size  0. 1% 

Clayton 

-0.5719 

-0.5643 

Frank 

-0.4925 

-0.4862 

FGM 

-0.6605 

-0.6712 

Gaussian 

-0.6005 

-0.5964 

Independent 

-1.5890 

-1.5824 

6.2  Sensitivities  of  Statistical  Moments 

Consider  the  same  2-D  highly  nonlinear  polynomial 
function  shown  in  Eq.  (31)  for  the  sensitivity  comparison  of 
statistical  moments.  Xj  and  X2  follow  Weibull  and  normal 
distributions,  respectively,  with  the  same  distribution 
parameters  as  shown  in  Table  7.  The  sensitivities  of  the  first 
two  statistical  moments,  mean  and  variance,  of  the 
performance  function  in  Eq.  (31)  are  obtained  using  Eqs.  (29) 
and  (30),  respectively.  Again,  1  million  MCS  samples  are  used 
for  the  sensitivity  calculation  using  Eqs.  (29)  and  (30),  that  is, 
K=  1,000, 000.  For  the  sensitivity  of  the  statistical  moment 
using  the  FDM,  the  statistical  moment  is  obtained  first  using 
the  MCS  with  250  million  samples.  Then,  by  perturbing  the 
mean  of  Xj  by  0.1%,  the  perturbed  statistical  moment  is 
obtained,  and  the  sensitivity  is  obtained  using  the  difference 
between  two  statistical  moments.  So,  for  a  problem  with  two 
random  variables  as  in  Eq.  (31),  the  sensitivity  analysis  using 
the  FDM  requires  three  MCSs  with  250  million  samples.  The 
computational  time  for  the  sensitivity  analysis  using  the  FDM 
is  750  times  that  of  the  proposed  sensitivity  analysis,  and  this 
difference  will  increase  as  the  number  of  random  variables 
increases. 

Table  12  compares  two  sensitivities  of  statistical  moments 
with  respect  to  jux  .  Regardless  of  the  copula  type  used, 
sensitivities  of  two  statistical  moments  obtained  using  the 
score  function  and  FDM  agree  very  well  with  each  other. 


Table  12.  Sensitivity  of  Statistical  Moments  with  respect  to  //, 


Mean 

Variance 

Copula 

Score 

Function 

FDM  with 
perturbation 
size  0.1% 

Score 

Function 

FDM  with 
perturbation 
size  0.1% 

Clayton 

11.8981 

11.6201 

-140.2574 

-136.1142 

Frank 

11.8093 

11.8681 

-109.1958 

-107.4361 

FGM 

11.2145 

11.1848 

-105.3351 

-104.3938 

Gaussian 

11.8488 

11.9547 

-122.4572 

-124.6348 

Independent 

11.2221 

11.2120 

-71.5755 

-71.1807 

6.3  Sampling-Based  RBDO  using  Proposed 
Stochastic  Sensitivity 

To  see  how  the  proposed  sensitivity  analysis  works  for  an 
RBDO  problem,  consider  a  2-D  mathematical  RBDO 
problem,  which  is  formulated  to 


minimize 
subject  to 


C(d)  =  - 


(dt  +d2  -10)2  (dj-^+10)2 


30 


120 


P(G,(X(d))  >  0)  <  P^  =  2.275%,  i  =  1  ~  3  (33) 


d"  <d<da. 


deR-andXeR2 


where  three  constraints  are  given  by 


G,(X)  =  1- 


XX 

20 


G,  (X)  =  -l  +  (7-6)2  +(7-6)3  -0.6x(7-6)4  +Z  (34) 
G3(X)  =  1-  80 


X2+8X2+5 


where 


0.9063  0.42261  |X] 

i  >  ,  which  are  drawn  in 
0.4226  - 0.9063 J  |Xj 


Figure  1.  The  properties  of  two  random  variables  are  shown  in 
Table  13,  and  they  are  assumed  to  be  correlated  with  the 
Clayton  copula  (x=0.5).  As  shown  in  Eq.  (33),  the  target 
probability  of  failure  ( PTFa' )  is  2.275%  for  all  constraints. 


Table  13.  Properties  of  Random  Variables 


Random 

Variables 

Distribution 

dL 

d° 

du 

Standard 

Deviation 

Xi 

Normal 

0.0 

5.0 

10.0 

0.3 

X2 

Normal 

0.0 

5.0 

10.0 

0.3 

Figure  1.  Shape  of  Constraint  and  Cost  Functions 


To  apply  the  proposed  sensitivity  analysis  to  an  RBDO 
problem,  accurate  surrogate  models  need  to  be  utilized.  For 
that  purpose,  the  D-Kriging  method  with  sequential  sampling 
proposed  by  Zhao  et  al.  [31]  is  used  because,  in  terms  of 
accuracy,  the  D-Kriging  method  outperforms  existing  methods 
such  as  the  polynomial  response  surface  method,  the  radial 
basis  function  (RBF)  method,  the  support  vector  regression 
(SVR)  method,  and  the  universal  Kriging  method.  Assuming 
the  constraint  functions  in  Eq.  (34)  are  not  given  analytically, 
which  is  usually  described  as  the  black-box  constraint,  the 
surrogate  models  for  the  constraint  functions  are  first 
generated  using  the  D-Kriging  method.  Then,  using  the 
generated  surrogate  models,  the  probabilities  of  failure  for 
three  constraints  and  their  sensitivities  are  obtained  using  the 
MCS  with  1  million  samples. 


Table  14  compares  the  numerical  results  of  five  different 
RBDO  methods.  The  first  three  results  are  obtained  from  the 
so-called  MPP-based  RBDO,  which  requires  the  sensitivity  of 
the  constraint  functions  for  the  MPP  search  and  design 
optimization.  This  MPP-based  RBDO  includes  the  FORM,  the 
DRM  with  three  quadrature  points,  which  is  denoted  as  DRM3 
in  Table  14,  and  the  DRM  with  five  quadrature  points,  which 
is  denoted  as  DRM5.  The  last  two  results  are  obtained  from 
the  sampling-based  RBDO,  which  uses  the  MCS  for  the 
estimation  of  the  probability  of  failure  and  its  sensitivity.  The 
sampling-based  RBDO  using  the  D-Kriging  method  is  the 
proposed  method,  and  for  the  comparison  of  the  accuracy  of 
the  proposed  method,  the  result  of  the  sampling-based  RBDO 
using  the  true  functions  given  in  Eq.  (34)  is  also  shown  in  the 
table. 


Table  14.  Comparison  of  Various  RBDOs  (PFar  =  2.275%  ) 


Optimum 

Design 

MCS 

Function 

Call 

Methods  Cost 

PFl ,% 

PFi  ,% 

MPP- 

FORM  -1.8742 

5.0026,  1.6165 

2.3022 

1.2835 

52+52 

Based 

DRM 3  -1.8794 

5.0315,  1.6050 

2.2912 

1.7496 

128+106 

RBDO 

DRM5  -1.8821 

5.0454,  1.5988 

2.2621 

2.0183 

146+102 

Sampling  D-Kriging  -1.8864 

5.0595,  1.5893 

2.2841 

2.3056 

57 

-Based 

RBDO 

Tn*e  -1.8853 
r  unction 

5.0541,  1.5918 

2.2912 

2.2791 

N.A. 

From  the  table,  it  can  be  seen  that  the  probability  of 
failure  of  the  second  constraint  (1.2835%)  at  the  optimum 
design  obtained  using  the  FORM  is  not  close  to  the  target 
probability  of  failure  (2.275%).  This  is  because  the  second 
constraint  is  highly  nonlinear  as  shown  in  Figure  1  and  the 
FORM  cannot  accurately  estimate  the  probability  of  failure  of 
highly  nonlinear  functions.  To  improve  the  accuracy  of  the 
probability  of  failure  at  the  optimum  design,  the  MPP-based 
DRM  with  three  or  five  quadrature  points  can  be  used;  Table 
14  shows  that  the  MPP-based  DRM  indeed  improves  the 
accuracy  of  the  probability  of  failure  at  the  optimum  design. 
Flowever,  to  obtain  a  more  accurate  optimum  design,  more 
quadrature  points  are  required,  such  as  the  DRM7,  etc.  The 
FORM  uses  52  function  evaluations  and  52  sensitivity 
calculations,  whereas  the  MPP-based  DRM  with  five 
quadrature  points  uses  146  function  evaluations  and  102 
sensitivity  calculations  to  obtain  the  optimum  design,  and  the 
number  of  function  evaluations  for  the  MPP-based  DRM  will 
be  increased  as  the  number  of  quadrature  points  increases. 

On  the  other  hand,  the  sampling-based  RBDO,  which  uses 
the  D-Kriging  method  and  the  proposed  stochastic  sensitivity 
analysis,  shows  very  accurate  optimum  design;  yet  it  requires 
only  57  samples  for  the  accurate  optimum  design.  Without  the 
sensitivity  of  the  performance  functions,  the  sampling-based 
RBDO  can  obtain  a  very  accurate  optimum  design  and  the 
optimum  design  is  very  close  to  the  optimum  design  obtained 
using  the  true  functions.  This  means  that  the  D-Kriging 
method  generates  very  accurate  surrogate  models  for  the  true 
functions.  From  this  example,  it  can  be  said  that  once  accurate 
surrogate  models  are  available,  the  sampling-based  RBDO 
using  the  proposed  sensitivity  analysis  yields  very  accurate 
optimum  designs  with  good  efficiency. 

More  detailed  discussion  on  the  sampling-based  RBDO, 
which  includes  higher  dimensional  engineering  applications 


such  as  M1A1  Abrams  tank  roadarm  [15],  needs  to  be  carried 
out  and  is  ongoing. 

7.  CONCLUSIONS 

The  stochastic  sensitivity  analysis  of  the  probabilistic 
constraints  and  statistical  moments  with  respect  to  mean 
values  of  correlated  random  variables  using  the  score  function 
and  MCS  is  carried  out  in  this  study.  Since  it  does  not  require 
the  sensitivity  of  performance  functions  or  even  the  sensitivity 
of  surrogate  models,  the  proposed  sensitivity  analysis  yields 
very  accurate  sensitivity  estimation  regardless  of  the  marginal 
and  copula  types  of  the  random  input  once  accurate  surrogate 
models  are  available.  Furthermore,  the  proposed  method  uses 
only  one  MCS  at  a  given  design  to  obtain  the  probability  of 
failure  and  its  sensitivity  or  statistical  moments  and  their 
sensitivities  simultaneously,  whereas  the  FDM  uses  N+l  MCS 
to  obtain  the  sensitivities  of  the  probabilistic  response,  where 
N  is  the  number  of  random  variables.  Thus,  the  proposed 
sensitivity  analysis  using  the  score  function  is  far  more 
efficient  than  the  FDM.  In  addition,  the  proposed  sensitivity 
analysis  is  more  accurate  than  the  FDM  when  the  same  MCS 
sample  size  is  used.  Hence,  the  proposed  stochastic  sensitivity 
analysis  combined  with  accurate  surrogate  models,  which  are 
obtained  in  this  paper  using  the  D-Kriging  method,  is 
recommended  for  RBDO  of  engineering  applications  where 
accurate  sensitivities  of  performance  functions  are  not 
available.  Numerical  examples  show  the  accuracy  of  the 
sensitivity  results  and  demonstrate  that  sampling-based  RBDO 
using  the  proposed  sensitivity  analysis  and  accurate  surrogate 
models  by  the  D-Kriging  method  yields  a  very  accurate 
optimum  design  with  good  efficiency. 
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